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Abstract. The linear and non-linear stability of sheared, relativistic planar jets is studied by means of linear 
stability analysis and numerical hydrodynamical simulations. Our results extend the previous Kelvin-Hemlholtz 
stability studies for relativistic, planar jets in the vortex sheet approximation performed by Perucho et al. (2004a, b) 
by including a shear layer between the jet and the external medium and more general perturbations. The models 
considered span a wide range of Lorentz factors (2.5 — 20) and internal energies (0.08 -60 c^) and are classified 
into three classes according to the main characteristics of their long-term, non-linear evolution. We observe a 
clear separation of these three groups in a relativistic Mach-number Lorentz-factor plane. Jets with a low Lorentz 
factor and small relativistic Mach number are disrupted after saturation. Those with a large Lorentz factor and 
large relativistic Mach number are the stablest, due to the appearance of short wavelength resonant modes which 
generate local mixing and heating in the shear layer around a fast, unmixed core, giving a plausible solution for 
the problem of the long-term stability of relativistic jets. A third group is present between them, including jets 
with intermediate values of Lorentz factor and relativistic Mach number, which are disrupted by a slow process 
of mixing favored by an efficient and continuous conversion of kinetic into internal energy. In the long term, all 
the models develop a distinct transversal structure (shear/transition layers) as a consequence of KH perturbation 
growth, depending on the class they belong to. The properties of these shear layers are analyzed in connection 
with the parameters of the original jet models. 
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1. Introduction 

■ Transversal structure in extragalactic jets could be the 
. natural consequence of current formation mechanisms 
' (see, e.g., Sol et al. 1989; Celotti & Blandford 2000), in 
which an ultrarelativistic presumably electron/positron 
outflow from the high latitude region close to the spinning 
black hole (and powered by, e.g., the extraction of energy 
from the hole) is surrounded by a mildly relativistic, elec- 
tron/proton, hydromagnetic outflow launched from the 
outer parts of the accretion disk. Recent numerical sim- 
ulations of jet formation from black hole magnetospheres 
(Koide et al. 1997) also lead to outflows with two-layered 
shell structure consisting of inner, fast gas pressure driven 
jets surrounded by slower, magnetically dominated winds. 
On larger scales, shear layers (with distinct kincmatical 
properties and magnetic fleld configurations) have been 
invoked in the past by several authors (Komissarov 1990, 
Laing 1996, Laing & Bridle 2002a, b) in order to account 
for a number of the observational characteristics of FR I 



Send offprint requests to: M. Perucho 



radio sources. The model of De Young (1993) to explain 
the FRI/FRII morphological dichotomy is based on decel- 
eration of the jet flow at the inner galactic core and the 
subsequent formation of turbulent shear layers in FRIs. 
Recently, Swain et al. (1998) found evidence of shear layers 
in FR II radio galaxies (3C353), and Attridge et al. (1999) 
have inferred a two-component structure in the parsec- 
scale jet of the source 1055-1-018. On the other hand, first 
simulations of radio emission from three-dimensional rel- 
ativistic jets with shear layers (Aloy et al. 2000) allowed 
several observational trends in parsec and kiloparsec jets 
to be interpreted: inhomogeneous distributions of appar- 
ent speeds within radio knots (Biretta et al. 1995); rails 
of low polarization along jets (as in 3C353; Swain et al. 
1998); top/down jet emission asymmetries in the blazar 
1055-1-018 (Attridge et al. 1999). Stawarz & Ostrowski 
(2002) have studied the contribution to the radiative jet 
output from turbulent shear layers in large-scale jets. 

Given all the pieces of theoretical and observational 
evidence concerning the ubiquity of shear layers in extra- 
galactic jets, it appears natural to analyze their influence 
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on the dynamics and stability of relativistic jets. An at- 
tempt to investigate the growth of the Kelvin-Helmholtz 
(KH) instability in some particular class of sheared, cylin- 
drical relativistic jets was pursued by Birkinshaw (1991). 
However, this study is limited to the ordinary and low- 
order reflection modes, and the domain of jet parame- 
ters considered involves only marginally relativistic flows 
(beam flow velocities < 0.1c; c is the speed of light) and 
non-relativistic (jet, ambient) sound speeds (< 0.01c). 
Other approaches to linear analysis of the stability of 
relativistic stratified jets (Hanasz & Sol 1996, Hardee & 
Hughes 2003) and sheared, ultrarelativistic jets (Urpin 
2002) have also been taken. Several recent works com- 
bine linear analysis and hydrodynamical simulations in 
the context of both relativistic jets (Rosen et al. 1999, 
Hardee 2000, 2001) and GRBs (Aloy et al. 2002). 



In this paper, we focus on the study of the evolution 
of relativistic (planar) fiows with shear layers through the 
linear and non-linear regimes relying on both linear, an- 
alytical stability analysis and hydrodynamical numerical 
simulations. The present work complements the one pre- 
sented in Perucho et al. (2004a, b; hereafter. Paper I and 
H, respectively), in which we characterized the effects of 
relativistic dynamics and thermodynamics in the devel- 
opment of KH instabilities in planar, relativistic jets in 
the vortex sheet approximation. We used a more general 
setup for simulations with the inclusion of a set of sym- 
metric (pinching) and antisymmetric (helical) sinusoidal 
perturbations in two dimensional slab jets and a thicker 
shear layer (~ 0.2 Rj) than that used to mimic vortex 
sheet evolution. The use of slab jets allows for inclusion 
of helical perturbations, which are known to be present 
in extragalactic jets. Moreover, two dimensional simula- 
tions provide the possibility of a much larger resolution 
than three dimensional ones. The aim of this work is to 
study the stability properties of jets depending on their 
initial parameters and the effect of shear layers in those 
properties. We used the temporal approach, which allows 
for larger resolution, and fixed two different grid sizes, de- 
pending on the thermodynamical properties of jets, which 
are neither directly related nor coupled to the wavelength 
of a specific mode, as was the case in Papers I and H. Jet 
parameters are based on those of previous papers for direct 
comparison. We generalized our results with simulations 
where only one antisymmetric mode is perturbed (similar 
to simulations in Paper I) , and simulations of cylindrical 
jets, where several modes are perturbed (as in simulations 
presented in this paper). 



The plan of this paper is as follows. In Sect. El we de- 
scribe the numerical simulations and present the parame- 
ters used in this paper. In Sect. 13 we describe our results 
concerning linear and nonlinear regimes of new simula- 
tions; we discuss them in Sect. 0] and present our conclu- 
sions in Sect. El 



2. Setup for numerical simulations 

The equations governing the evolution of a relativistic 
perfect-fluid jet (see Paper I) are 



2f , P 



^ + (..V). 



^ + V.(7.) 



= 0, 



djpPo ^) 
dt 



0. 



(1) 



(2) 



(3) 



In the preceding equations, c is the speed of light, po the 
particle rest mass density, p stands for the relativistic den- 
sity which is related to the particle rest mass density and 
the specific internal energy e, hy p — pq{1 + e/c^). The 
enthalpy is defined a,s w = p + p/c^, and the sound speed 
is given by = Tp/w, where F is the adiabatic index. The 
relation between pressure and the specific internal energy 
is p = (r — l)epo- The velocity of the fiuid is represented 
by V, and 7 is the corresponding Lorentz factor. The oper- 
ator d /dt appearing in Eq. Q is the standard Lagrangian 
time derivative. 

The steady initial flow is 2D-planar and symmetric 
with respect to the a; = plane. The flow moves in the pos- 
itive z direction. Our simulations were performed in the 
so-called temporal approach, as in Papers I and II. In this 
approach, the evolution of perturbations in a peridiocal 
slice of an inflnite jet is followed along the time. In order 
to study the eflects of shearing in the development of the 
instability, we assumed a continuous transition between 
the jet and the ambient (the same as the one considered 
by Ferrari et al. 1982). The profiles of the axial velocity 
and proper rest-mass density across this transition layer, 
Vz{x) and po{x), respectively, are given by 



cosh{x /RjY 



P0(X) = po,a 



POm — PQ,j 

coshix/RjY' 



(4) 



(5) 



In the previous expressions, represents the fiuid flow 
velocity in the jet axis, whereas poj and po,a are the 
proper rest-mass density at the jet axis and in the ambi- 
ent, respectively. The exponent m controls the shear layer 
steepness; in the limit m 00 the configuration tends to 
the vortex-sheet case. In our present calculations we have 
used m = 25, corresponding to a shear layer of thickness 
0.17 about twice that used in Papers I and II, in order 
to mimic the vortex sheet limit. From now on all quanti- 
ties representing the jet will be assigned the 'j' subscript 
and the quantities representing the ambient medium will 
be assigned the 'a'. 

Following conclusions given in the Appendices of both 
Paper I and Paper II, the numerical resolution used was 
256cells/i?j in the transversal direction times 32cells/i?j 
in the axial direction. Note that we reduced the transversal 
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Model 


7 


Ej 


A2.5 


2.5 


0.08 


B2.5 


2.5 


0.42 


D2.5 


2.5 


60.0 


B05 


5 


0.42 


D05 


5 


60.0 


AlO 


10 


0.08 


BIO 


10 


0.42 


DIO 


10 


60.0 


B20 


20 


0.42 


D20 


20 


60.0 



0.008 
0.042 
6.000 
0.042 
6.000 
0.008 
0.042 
6.000 
0.042 
6.000 



0.18 
0.35 
0.57 
0.35 
0.57 
0.18 
0.35 
0.57 
0.35 
0.57 



0.059 
0.133 
0.544 
0.133 
0.544 
0.059 
0.133 
0.544 
0.133 
0.544 



P 

0.0027 
0.014 
2.000 
0.014 
2.000 
0.0027 
0.014 
2.000 
0.014 
2.000 



0.11 0.11 
0.14 0.15 
0.87 0.90 
0.14 0.15 
0.87 0.90 
0.11 0.11 
0.14 0.15 
0.87 0.90 
0.14 0.15 
0.87 0.90 



12.5 
6.12 
3.29 
13.2 
7.01 
54.2 
26.9 
14.2 
54.0 
28.5 



ko ki 



0.39 
0.39 
0.78 
0.39 
0.78 
0.39 
0.39 
0.78 
0.39 
0.78 



0.78 
0.78 
1.57 
0.78 
1.57 
0.78 
0.78 
1.57 
0.78 
1.57 



fc2 ks 

1.18 1.57 

1.18 1.57 

2.36 3.14 

1.18 1.57 

2.36 3.14 

1.18 1.57 

1.18 1.57 

2.36 3.14 

1.18 1.57 

2.36 3.14 



Table 1. Equilibrium parameters of different simulated jet models. The meaning of the symbols is as follows. 7: jet 
flow Lorentz factor; e: specific internal energy; Cg'. sound speed; p: pressure; p: jet-to-ambient relativistic rest mass 
density contrast; 77: jet-to-ambient enthalpy contrast; Mji jet relativistic Mach number; /co,i,2,3: excited longitudinal 
wavenumbers. Labels a and j refer to ambient medium and jet, respectively. All the quantities in the table are expressed 
in units of the ambient density poa, the speed of light c, and the jet radius Rj. 



resolution with respect to the simulations in Papers I and 
II. One reason for that were computational time limita- 
tions, as now our grids are twice as large in the transversal 
direction as those used in Paper I, since we are now com- 
bining symmetric and antisymmetric structures. However, 
in the present simulations, transversal resolution is not 
as critical as in the previous works, since we are not in- 
terested in mimicking the evolution of instability in the 
vortex sheet limit and therefore do not have steep shear 
layers. The saving in transversal resolution allowed us to 
double axial resolution, which affects the non-linear results 
(see Appendix in Paper II) . The physical sizes of grids are 
8 Rj axially times 6 Rj transversally for hot jets (models 
D, see Table^ and 16 Rj axially times 6 Rj transversally 
for cold jets (models A and B in Table ^1. The different 
axial size is due to hot models having shorter unstable 
modes; see Sect. 13.11 where we show linear problem solu- 
tions for one cold and one hot jet. 

Previous to these simulations, we performed several 
representative runs (BOS, DOS, B20, D20) with the aim of 
studying the evolution of models under single eigenmode 
perturbations in planar antisymmetric geometry, exciting 
the first reflection antisymmetric mode at its peak growth 
rate. These simulations were used to check the consistency 
of the numerical results in the linear phase with the linear 
stability analysis for relativistic, sheared flows. Discussion 
of the evolution of these models through the non-linear 
regime can be found in Appendix A. Tests were also per- 
formed in order to assess the difference in the evolution 
of linear and non-linear regimes using a general sinusoidal 
perturbation (to be used in this paper) and a superpo- 
sition of eigenmodes, as done with the first body mode 
alone in Paper I. Results showed that structures and qual- 
itative properties of the resulting flow were basically the 
same. This fact confirms that general perturbations excite 
eigenmodes of the system. 

The parameters used in the simulations arc shown in 
Tabled We swept a wide range in Lorentz factors (from 
2.S to 20) and in specific internal jet energies (from 0.08c^ 



to 60c^) in order to obtain a global view of the response of 
different initial data sets to perturbations. All the models 
correspond to a single-component ideal gas with adiabatic 
exponent F = 4/3. These parameters were chosen in order 
to study the stability regions found in Paper II: Class I 
for cold and slow models, which were deeply mixed and 
mass loaded; Class II for hot and fast jets, which were 
slowly mixed in the nonlinear regime, progressively losing 
their axial momentum; Class III for hot and slow jets, 
with properties between Classes I and II, and Class IV for 
cold/warm and fast models, which were the stablest in the 
nonlinear regime. We performed simulations for models 
B05\ BIO, B20, DOS, DIO, and D20 of Paper II, and 
added A2.5 (same thcrmodynamical properties as AOS in 
Paper II), AlO, B2.5, and D2.5. Models A2.5, B2.5, 
and BOS correspond to regions of Class I jets. Models 
DIO and D20 correspond to Class II, D2.S and DOS to 
Class III, and A10,B10 and B20 belong to Class IV. 

Perturbations we applied adding the following sinu- 
soidal form to transversal velocity, Vx{x, z): 

Vx = 

( X! sin((n+ l)knZ + v?„)sin^((n+ l)7rx)— j + 
-j^ I X! '5in((m -I- 1) z + (p,n) sin2((TO + 1) n x) \ , (6) 

\m=0 ) 

where Vx\{^ lO^**) is the amplitude given to the pertur- 
bation, fcfn.n are the wavenumbers of the grid box (so that 
(n -t- 1) fc„ and (m -f 1) fcm stand for the harmonics of the 
symmetric (pinching) and antisymmetric (helical) modes, 
respectively), and and are random phases given to 
each mode. In our simulations, four symmetric {N = 4) 
plus four antisymmetric modes (M = 4) were excited, i.e., 
the fundamental mode of the box and the first three har- 
monics. 

^ Boldface will be used for new simulations in order to dif- 
ferentiate them from those in Paper II with the same name. 
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Numerical siniulations were performed using a finite- 
difference code based on a high-resolution shock-capturing 
scheme which solves the equations of relativistic hydrody- 
namics written in conservation form (Marti et al. 1997 
and references therein). Before performing the simula- 
tions, several improvements were made in the numerical 
code. In particular, the relativistic PPM routines (Marti 
& Miiller 1996) were properly symmetrized. The code was 
recently parallelized using OMP directives. Simulations 
were performed with 4 processors on SGI 2000 and SGI 
Altix 3000 machines. 

3. Results 

3.1. Linear regime 

3.1.1. Perturbation theory 

We introduced an adiabatic perturbation of the form 
oc g{x) ex.p(i{kzZ — uit)) to the flow equations (|1I3|I . u 
and kz (kx) being the frequency and wavenumber of the 
perturbation along (across) the jet flow. We followed the 
temporal approach, in which perturbations grow in time, 
with real wavenumbers and complex frequencies, with the 
imaginary part being the growth rate. By linearizing the 
equations and eliminating the perturbations of rest mass 
density and flow velocity, a second-order ordinary differen- 
tial equation for the pressure perturbation Pi , is obtained 
(Birkinshaw 1991, Aloy et al. 2002) 



A" 



Pefi 



^ - vozkz 

{UJ - VQzkz f 



Pe,0 + Pq 

{kz - ujvoz 



(7) 



Pi 



where pe.o is the energy-density of the unperturbed model, 
Pq the pressure, vqz the three-velocity component, 70 — 
Xj \J\ — the Lorentz factor, and c^.o the relativistic 
sound speed. The prime denotes the x-derivative. Unlike 
the vortex sheet case, in the case of a continuous velocity 
profile, a dispersion relation cannot be written explicitly. 
The equation ((21 is integrated from the jet axis, where 
boundary conditions on the amplitude of pressure pertur- 
bation and its first derivative are imposed 



P^[x = 0) = 1, Pi'(a; = 0) = (sym. modes), 
Pi[x = 0) = 0, P[{x = 0) = 1 (antisym. modes). 



(8) 



Solutions satisfying the Sommerfcld radiation conditions 
(no incoming waves from infinity and wave amplitudes de- 
caying towards infinity) are found with the aid of the nu- 
merical method, based on the shooting method (Press et 
al. 1997) proposed in Roy Choudhury & Lovelace (1984). 

Linear stability analysis was performed for all models 
presented in Sect. |21 in the symmetric and antisymmet- 
ric cases. Figure ^ shows examples of solutions for the 
linear problem for sheared jets in models BOS and D20. 
Top panels in Fig. ^ show the (real part of) frequency as 
a function of wavenumber, and bottom panels show the 



imaginary part of frequency or growth rate, defined as the 
inverse of the time needed by a given mode to e-fold its 
amplitude, also as a function of wavenumber. Each mode 
is defined by its wavelength, frequency and growth rate. 






Fig. 1. Solutions for the linear perturbation differential 
equation. Left panel: Antisymmetric solution for model 
BOS. Right panel: Symmetric solution for model D20. 
Vertical lines stand for the perturbed wavenumbers in the 
numerical simulations (from left to right: k^, ki, ^2 and 
fca). Let us note that the fundamental mode does not ap- 
pear in the right panel, as its growth rates are lower than 
the scale of the plot. 

From these results, we note that the individual reflec- 
tion mode solutions of the shear problem present lower 
growth rates for most wavenumbers, especially in the large 
wavenumber limit, than do the corresponding solutions 
in the vortex sheet case. This behaviour was noticed for 
the first time by Ferrari et al. (1982) for the first and 
second reflection modes in the non- relativistic limit. The 
growth rate curves corresponding to a single rix-th re- 
flection mode consists of a broad maximum at higher 
wavenumbers and a local peak which is placed in the 
low wavenumber limit, near the marginal stability point 
of a chosen reflection mode. Regarding the relativistic 
case, while in the vortex-sheet limit the small wavenum- 
ber peaks for individual modes are relatively unimportant 
(since the maximum growth rates at these peaks are lower 
than those of other unstable modes), while in the presence 
of the shear-layer they display high growth rates for high 
order body modes. Therefore we shall call these peaks the 
shear layer resonances. In Fig. [21 we show the solution for 
four specific modes of model D20, from Fig.^ Low order 
body modes do not show high peaks at maximum unstable 
wavelengths, whereas high order body modes show peaks 
(shear layer resonances) at this maximum wavelength and 
do not present broad maxima. 

From Eq. Q we see that radial structure of pertur- 
bations depends on physical parameters of the flow, as 
well as on the given frequency and axial wavenumber 
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Fig. 2. Specific modes of solutions shown in Fig.^ sym- 
metric solution for model D20. Dotted line: first body 
mode; dashed: second body mode; dash-dot: twentieth 
body mode; dash-triple dot: twenty-fifth body mode. 
Arows point to both the broad maxima and the small 
wavenumber peaks present in every single mode. Low 
wavenumber peaks of high order body modes show higher 
growth rates and are thus defined as (shear layer) reso- 
nances. 



of a given perturbation. Resonances are determined by 
this transversal structure, and therefore we should expect 
changes in their characteristics depending on the proper- 
ties of the shear layer and physical parameters: i) a de- 
crease of the jet Lorentz factor reduces the dominance of 
the growth rates of resonant modes with respect to or- 
dinary and low order reflection modes; ii) a decrease in 
the specific internal energy of the jet causes resonances to 
appear at longer wavelengths; iii) further widening of the 
shear layer reduces the growth rates and the dominance 
of the shear-layer resonances, suggesting that there is an 
optimal width of the shear layer that maximizes the ef- 
fect for a given set of jet parameters; the largest growth 
rate of resonant modes moves to smaller wavenumbers 
and lower order reflection modes; iv) perturbations with 
wavenumber higher than some limiting value (that de- 
creases with the shear layer width) are significantly dimin- 
ished (short- wavelength cut-off), consistent with previous 



non-relativistic results (Ferrari et al. 1982). The discov- 
ery of the shear layer resonances and their potential role 
in the long-term stability of relativistic jets is reported in 
Perucho et al. (2005). 

3.1.2. Simulations 

Table 12 summarizes the properties of the linear phase in 
our simulations. The left part of the Table (colums 2-9) 
gives the theoretical growth rates of the perturbed wave- 
lengths, taken at the vertical lines in Fig. ^ The last col- 
umn gives the values of the growth rate corresponding to 
the dominant wavelength as deduced from Fourier anal- 
ysis of the transversal profiles of the rest mass density 
distribution in the jet. Note, however, that Fourier anal- 
ysis can only give us information about wavelengths, but 
cannot distinguish between symmetric and antisymmet- 
ric modes. The growths of pressure, axial, and transversal 
velocity perturbations along the simulations are shown in 
Fig. El 

Comparison of the evolution through the linear phase 
of the different models in the numerical simulations and 
from the linear stability analysis is summarized as follows: 

— A2.5: modes with longer wavelengths are faster grow- 
ing, and their Fourier amplitudes are consistently 
larger than those for shorter modes in the simulation. 
Growth rate found in the simulation is close to the one 
expected from linear stability analysis. 

— B2.5: first (fci) and second (fc2) harmonics of the box 
have larger amplitudes in the Fourier analysis and, 
therefore, dominate the linear regime. Linear stabil- 
ity analysis gives fci, k2 and as the fastest growing 
modes, with a similar growth rate to that found in 
the simulation. However, modes have smaller am- 
plitudes. 

— D2.5: found growth rate for the simulation is close 
to that of fci and fc2 modes, which is confirmed by 
Fourier analysis. Antisymmetric fcs mode might grow 
with slower rates than theory predicts due to numerical 
viscosity that affects shorter modes more than longer 
ones. 

— BOS: Fourier analysis shows competition between fun- 
damental and first harmonics of the box (fco and fci, 
respectively). This, as well as the value of the mean 
growth rate, is confirmed by the linear stability anal- 
ysis. The second harmonic of the box (fc2) is damped. 

— DOS: according to Fourier analysis, fci and fc2 modes 
dominate evolution in the linear regime. The growth 
rate is close to that of the symmetric fci mode, de- 
spite the fact that symmetric fc2 and antisymmetric 
fcs present faster growth rates, so they appear to be 
damped. 

— AlO: Fourier analysis shows that longer modes dom- 
inate, in agreement with linear stability analysis. 
However, the growth rate in the numerical simulation 
is two times smaller than predicted. We also observe in 
Fourier analysis that very short modes, excited as har- 
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Model 




fco 


w. 


,ki 


Wi 




w. 


,k3 


Dominant 


Wi 




Symm. 


Antis. 


Symm. 


Antis. 


Symm. 


Antis. 


Symm. 


Antis. 






A2.5 


0.036 


0.032 


0.038 


0.037 


0.034 


0.036 


0.031 


0.034 


ko 


0.030 


B2.5 


0.042 


0.056 


0.070 


0.052 


0.066 


0.084 


0.073 


0.080 


ki,k2 


0.070 


D2.5 


0.046 


0.160 


0.131 


0.182 


0.210 


0.194 


0.142 


0.256 


k2,ki 


0.200 


BOS 


0.037 


0.035 


0.037 


0.044 


0.036 


0.038 


0.034 


0.035 




0.035 


D05 


0.068 


0.063 


0.085 


0.063 


0.100 


0.068 


0.068 


0.110 


ki,k2,ko 


0.080 


AlO 


0.009 


0.009 


0.006 


0.006 


0.005 


0.006 


0.006 


0.006 


ko* 


0.004 (0.005) 


BIO 


0.022 


0.018 


0.019 


0.021 


0.018 


0.017 


0.013 


0.013 


ko 


0.020 


DIG 


0.034 


0.038 


0.041 


0.037 


0.044 


0.034 


0.051 


0.035 


ki,k2 


0.040 


B20 


0.011 


0.010 


0.009 


0.010 


0.007 


0.007 


0.009 


0.010 


ko* 


0.006 (0.008) 


D20 


0.018 


0.018 


0.020 


0.017 


0.022 


0.017 


0.027 


0.028 


ki,ko 


0.016 



Table 2. Dominant modes in the linear phase of the numerical simulations. Wi^kj- maximum growth rate for the 
j-th overtone wavenumber excited in the simulation (see Table ^ derived from linear stability analysis. Left columns: 
symmetric mode; right columns: antisymmetric one. The dominant mode refers to the mode with the largest amplitude 
in rest mass density perturbation as derived from Fourier analysis of the box, written from larger to smaller amplitude 
when more than one is present, wf. fitted pressure perturbation growth rate for the linear regime in the simulation. 
Growth rate values are in c/Rj units. *: Models where irregular growth affects the evolution (see text). 



monies of perturbed wavelengths, become important 
by the end of the linear regime. 

— BIO: fco modes dominate, as predicted by linear sta- 
bility analysis. 

— DIO: as in models D2.5 and DOS, fci and fc2 have 
larger amplitudes in Fourier analysis, but the smaller 
wavelength modes (fcs) are damped with respect to the 
predictions of linear stability analysis. 

— B20: longer modes dominate the linear evolution, in 
agreement with linear analysis, but the growth rate in 
the numerical simulation is 1.5 times smaller than pre- 
dicted. After some time, short, fast modes, like those 
appearing in model AlO, become dominant and lead 
to a smooth transition to the non-linear regime. 

— D20: long modes present larger amplitudes with pre- 
dicted growth rates up to the moment when shorter 
modes reach larger amplitudes, the same effect as 
found in models AlO and B20. 

It is observed in several simulations (e.g., B2.5, BOS, 
D2.S, DOS, DIO) that modes with similar or even slightly 
higher growth rates than those dominating in simulations 
present smaller amplitudes in the linear regime. It hap- 
pens usually for shorter modes (typically k2, fca), so it 
may be caused by numerical viscosity, for less cells are in- 
volved in one wavelength. However, the way in which we 
perturb the jet may also favor the dominating growth of 
certain modes starting with a larger amplitude. We added 
a general sinusoidal perturbation, so the input amplitude 
of the perturbation at a given wavelength is shared in a 
random way among all the modes present at that wave- 
length. This makes some modes start their growths with 
smaller amplitudes, as we could see in the Fourier analysis 
of different models. Initial low amplitudes are more prob- 
able for short wavelength modes, as more eigenmodes are 
present at a given wavenumber in this range (see Fig.^. 
From an initial lower amplitude, and taking into account 
that they have similar growth rates to other modes, they 



Model 


''^i .max 






B05 


0.052 






D05 


0.11 






AlO 


0.013 


0.017 


0.009 


BIO 


0.035 






DIO 


0.057 






B20 


0.026 


0.036 


0.036 


D20 


0.035 


0.070 


0.047 



Table 3. Growth of resonant modes. Models which 
present a global maximum growth rate (according to the 
linear analysis) for all resonant modes (i.e., at any wave- 
length) above the growth rates of the perturbed modes 
are listed. Wi^max- maximum growth rate for all resonant 
modes from linear analysis; Wi^p^^^ : fitted growth rates of 
pressure and perpendicular velocity perturbations for the 
fast growth linear regime in the simulation, only for those 
simulations where it occurs; Wi \\: same as Wi^p^y^^ for axial 
velocity. All values are in c/Rj units. 



grow with smaller amplitudes for the rest of the linear 
phase. 

Models AlO and B20, marked with an asterisk in 
Table |2 have fitted growth rates in the first part of 
the linear regime below the predicted values. Note that 
these models have the lower growth rates. After this ini- 
tial phase, short harmonics start dominating the linear 
growth. 

We have observed the appearance of fast growing, very 
short modes in models AlO, B20, and D20, which are 
clearly associated to the resonant modes presented above 
in the previous section and which could have been ex- 
cited as harmonics of the initially perturbed wavelengths. 
The same kind of resonant mode might have developed 
in model C20 of Paper II and caused the irregular linear 
growth found with respect to the rest of models. These 
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Fig. 3. Evolution of the relative amplitudes of perturbations. Dotted line: pressure perturbation ((p„ 



Po)/po)- 



Dashed hue: longitudinal velocity perturbation in the jet reference frame (0.5 (u^ ^^^ — v'^^^^). Dash-dotted hne: 



perpendicular velocity perturbation in the jet reference frame (0.5 (v'^ 
horizontal axis. The search for maximum (pmax , v'^ , v'^ 
to those numerical zones with jet mass fraction larger than 0.5. 



vi ^„^, vi ^„^) and minimum (v'^ 



J). Note the different scales in the 
i,mm> '"'z,min) valucs wcrc restricted 



modes generate a rich internal structure in the jet due 
to their large perpendicular wavenumber or, equivalently, 
short perpendicular wavelengths (characteristic of high or- 
der body modes). A direct comparison between the struc- 



ture generated by these resonant modes in the numerical 
simulations and that coming from linear stability anal- 
ysis can be seen in Fig. ^ In this figure, we display one 
snapshot from model D20 and the theoretical counterpart 
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Fig. 4. Upper panels: pressure (left) and perpendicular velocity perturbation (right) at late stages of linear phase 
(model D20). Lower panels: pressure (left) and perpendicular velocity perturbation (right) corresponding to one 
resonant mode from linear analysis. The linear, grey scale is arbitrary. Amplitudes are maxima at the shear layer, 
hence the name of shear layers resonances given to these modes. Oblique waves in upper panels are the result of longer 
wavelength perturbations, not present in the bottom panels. 



using one of those resonant modes. The upper plots cor- 
respond to the numerical simulation in the linear regime, 
where the signature of the initial perturbations (fco,fci,fc2 
and fca) are the oblique waves seen outside the jet. As 
seen in these plots, resonant modes grow to amplitudes 
larger than those of the long waves, as indicated by the 
black/white scale saturating precisely on the interface. 
The lower plots represent the theoretical structure that 
we would find if we had only excited a resonant mode and 
that is fairly similar to the one appearing in the shear 
layer of the simulated jet. Let us point out, however, that 
it is difficult to identify the exact mode in the simulation, 
as the resonant modes overlap so much (see Fig. 1) and 
it may happen that what we see is the structure resulting 
from a combination of competing resonant modes. 

According to the linear stability analysis, resonant 
modes have the highest growth rates in high Lorentz factor 
jets and, among them, in colder jets. This could be the rea- 
son why they only appear in simulations of models AlO, 
B20, and D20. Table El collects the models which present 
maximum growth rate for all resonant modes (i.e., at any 
wavelength), found in the solutions to the linear problem, 
above the growth rates of the perturbed modes. Maximum 
growth rates for resonant modes in those models where 
they have been found, along with the fitted growth rates 
in the simulation, are listed. Typically, the growth rates 



from the numerical simulations are about 1.4 — 2.0 times 
higher. This difference remains unexplained, but it could 
be caused by second-order effects, like interaction between 
modes. 

Summarizing, two kinds of linear growth are found in 
these simulations, one dominated by longer modes typi- 
cal of slower jets and another one where short, fast modes 
appear. This difference is important, for the transition to 
the non-linear evolution depends critically on the domi- 
nant modes at the end of the linear regime. 

Table 01 shows the times at which linear phase ends. 
As the end of the linear regime we selected the moment 
when one of the variables (usually axial velocity) changes 
its slope (departs from the linear growth, see definitions in 
Paper I). On the other hand, we noticed that the longitudi- 
nal velocity perturbation grows linearly up to values close 
to the speed of light and then beyond the sound speed. 
This means that shocks should form at the end of linear 
phase, as it is the case; see Fig. 3 in Paper I, where we ob- 
serve weak shocks starting to appear as conical structures. 
We could have selected the end of the linear regime as the 
moment when these shocks start to appear. This would 
relate the end of the linear regime directly to the internal 
sound speed. We see that colder jets have longer linear 
phases than hot ones, due to smaller typical growth rates 
in the former, tun times are larger than those in Paper I, 
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Model 


ilin 


^mcx 


^mix 


^sat 


^pcak 


^pcak 


A2.5 


225 


250 


300 


340 


340 


100 


B2.5 


110 


125 


140 


150 


165 


20 


D2.5 


40 


45 


50 


50 


50 


2 


BOS 


220 


275 


300 


280 


330 


70 


D05 


105 


125 


110 


130 


140 


2 


AlO 


725 


— 


— 


— 


— 


3 


BIO 


400 


520 


500 


500 


540 


100 


DIG 


205 


260 


220 


290 


300 


4 


B20 


475 




520 


550 


590 


4 


D20 


275 


510 


300 


275 


320 


2 



Table 4. Times for the difFercnt phases in the evolution 
of the perturbed jet models, tn^: end of Unear phase (the 
amphtudes of the different quantities are not constant any 
longer). tsa.t- end of saturation phase (the amplitude of the 
transverse speed perturbation reaches its maximum), iniix^ 
the tracer starts to spread, tpcak'- the peak in the amplitude 
of the pressure perturbation is reached, tmc^: the jet has 
transferred to the ambient 1% of its initial momentum. 
Apoak^ relative value of pressure oscillation amplitude at 
the peak of pressure perturbation (see Fig.|31). 



as growth rates are reduced by the presence of the shear 
layer. Model AlO presents the longest linear phase. 

3.2. Saturation and transition to nonlinear regime 

Saturation of perturbations is reached (see Paper I) when 
perpendicular velocity cannot grow further in the jet refer- 
ence frame due to the speed of light limit. Saturation times 
tsat for the different models are listed in Table 01 In this 
phase, structures generated by dominating modes become 
visible in the deformations of the jet. In Fig. [S] we show 
snapshots of three models (B2.5, DOS, and DIO) at satu- 
ration time where mode competition derived from Fourier 
analysis is clearly observed. Asymmetric structures appear 
as a result of several symmetric and antisymmetric modes 
with large amplitudes. 

In Paper II we also discussed how at the end of the 
saturation phase nearly all the simulations lead to a sharp 
peak in the pressure oscillation amplitude. These peaks 
are also seen in the present simulations (see Fig. The 
relative values of pressure oscillation amplitude at the 
peak Apeak and the corresponding times tpoak are listed 
in Table 01 The values of Apeak were connected with the 
non-linear evolution of the flow. Those cases in which 
Apeak > 70 developed a shock in the jet/ambient interface 
followed by the sudden disruption of the jet. From TableQ] 
we see that peak values in the present simulations are in 
general qualitatively the same as the corresponding ones 
in Paper I. Colder and slower jets have larger peaks and 
hence suffer stronger shocks after saturation. The main 
difference between the values in this paper and those pre- 
sented in Paper I appears for models B20 and D20, where 
shock strength is much smaller due to the appearance of 
resonant, stabilizing modes, as we discuss next. 



The parallel and perpendicular wavelengths of the 
shear-layer resonant modes, Az and Xx, respectively, are 
both small (Az < Rj) with A^; <C Az. Therefore their 
wavevectors are almost perpendicular to the jet axis so 
the waves propagate from the shear layer towards the jet 
axis. On the other hand, the resonant modes have high 
growth rates, exceeding the growth rate of other modes, 
so they start to dominate in the evolution. Subsequently, 
the resonant modes saturate as soon as the flow velocity 
oscillation amplitude approaches the speed of light. As the 
maximum amplitude is reached, the sound waves steepen 
while travelling towards the jet axis and form shock fronts 
on the leading edges of wave profiles Dissipation of the os- 
cillation energy of resonant modes in shocks changes the 
background flow, so that the amplification conditions of 
the longer wavelength modes change during the course of 
time, reducing the value of Apeak and preventing the for- 
mation of a strong shock. 

Finally, as found in Paper II, the generation of the 
shock wave at the jet/ambient interface is imprinted in 
the evolution of the maxima of the transversal Mach 
number of the jet with respect to the unperturbed am- 
bient medium. This quantity is defined as Mj^j_ = 
lj,±Vj,±/{jcsa'^sa), with 7j^_L and 7c^„ the Lorentz factors 
associated to Vj^± and Csa, respectively. A value signifi- 
cantly larger than 1 around tpeak points towards a super- 
sonic expansion of the jet at the end of the saturation 
phase. This magnitude is shown in Fig. |B1 We observe a 
clear inverse tendency of the peak value of this magnitude 
from colder to hotter and from slower to faster jets, with 
the exception of AlO with respect to BIO and DIO, due 
to the presence of the resonant stabilizing modes prevent- 
ing the formation of a strong shock. 

It is important to note that models with Apeak > 10 
(A2.5, B2.5, B05, and BIO) coincide with those develop- 
ing larger transversal Mach numbers, see top panels of Fig. 
[Sjfor model B2.5, where pressure maxima are observed at 
the jet center and in the interaction of the growing wave 
with the ambient. 

3.3. Fully non-linear regime 

In Paper II, the non-linear evolution of the instability in 
the different models was characterized through the pro- 
cesses of jet/ambient mixing and momentum transfer. In 
Fig. 13 we show the width of the mixing layer as a func- 
tion of time for all the models. The times at which mixing 
starts in the different models ^mix appear listed in Tabled 
In all cases these times are around tgat- Generally, mod- 
els developing wider shear layers are also more mixed; i.e., 
the amount of mass in zones with jet mass fraction strictly 
different from and 1, is higher. We observe that those 
models with larger values of Apeak (lower Lorentz factor 
and colder jets) develop wider layers (> 5i?j) soon after 
saturation due to turbulent mixing induced by the shock, 
while those models where resonant modes appear do not 
show strong mixing with the ambient. Models BIO and 



10 



M. Perucho et al.: Nonlinear stability of relativistic sheared planar jets 

Time= 150R/c 



•t.3SB-01 



Log Pressure 





2 4 6 3 W 12 14 16 



2 4 e B 10 12 14 16 



Time = 130 R/C 

Log Pr&ssufe 



1.036+01 



1.77B+00 



-4 

3.03S-01U -6 




7.196+00 



4.1Oa+00 



1.00S+00 



S * 6 9 10 12 1* ie 

Time = 290 R/c 

Log Pressure 



2.7Se-01 



B.45a-a3 





6 




4 - 




2 -r 




- 






1 


t 







I.Bis+OI 



B.S4e+00 



LOOs+aa 



a 2 4 e A 10 12 u 16 



Flaw Lofentz Factor 




O 2 4 S B 10 12 14 IS 



Flow Lorentz Factor 




2 4 e B 10 IS^ 14 16 



Fig. 5. Snapshots of logarithm of pressure (left) and Lorentz factor (right) for models B2.5 (upper panels), DOS 
(center panels) and DIO (bottom panels) at tsat, where irregular structures caused by mode competition are observed. 



DIO undergo a mixing process, though slower than the 
former. 

Figure IHl shows the fraction of initial axial momentum 
kept by the jet as a function of time. Axial momentum is 
lost first through sound waves forming the linear pertur- 
bations and second, but more important, through shocks 
themselves and by subsequent mixing, which implies shar- 
ing of momentum with the ambient medium. Correlation 
with Fig. d is remarkable. Models developing wide mix- 
ing layers coincide with those losing more than 50% of 
their initial axial momentum just after saturation; models 
BIO and DIO share their momentum with the ambient 
medium continuously in the non-linear regime; and mod- 
els where resonant modes dominate saturation keep almost 
all their initial momentum by the end of the simulations. 
Results derived from Fig. |H| are corroborated by Fig. In 



the latter we plot the total transversal momentum in the 
jets normalized to the corresponding longitudinal momen- 
tum. Transversal momentum in the jet (initially zero) is 
generated through turbulent motions and continuous con- 
version of kinetic into internal energy. The value of the 
normalized transversal momentum at a given time is an 
indication of how far from equilibrium the jet is. We ob- 
serve that colder and lower Lorentz factor models present 
strong peaks at tgat , coincident with the triggering of the 
shock and the sudden transfer of longitudinal momentum 
seen in the previous plot: Those models where resonant 
modes appear barely generate any transversal momen- 
tum, and models BIO and DIO do not present strong 
peaks at saturation but display a steady transmission of 
the transversal momentum through the non-linear regime 
(see Fig. inj, implying continuous loss of energy. 
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Fig. 6. Transversal Mach number in simulations (see text for definition). Solid line: 7 = 2.5; dotted line: 7 = 5.0; 
dashed line: 7 = 10.0; dash-dot line: 7 = 20.0. The thinnest line is for model A and thickest for model D, with B 
in between, the left panel shows models A2.5, B2.5, BOS, and BIO, and the right panel shows models AlO, B20, 
D2.5, D05, DIO, and D20. 
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Fig. 7. Evolution of the mean width of the jct/ambicnt 
mixing layer (for tracer values between 0.05 and 0.95) with 
time for all simulations. Lines represent the same models 
as in Fig. El A value of 5 Rj for the width of the mixing 
layer (horizontal dashed line) serves to classify the evolu- 
tion of the different models. 



Fig. 8. Evolution of the normalized total longitudinal mo- 
mentum in the jet as a function of time. Lines represent 
the same models as in Fig. El The long-dashed horizontal 
line identifies those models transferring more than 50 % 
of the initial jet momentum to the ambient. 



Panels showing several physical quantities for all mod- 
els at the end of simulations are presented in Figs. I10I19I 
Colder and slower models (A2.5, B2.5, and BOS) show 
turbulent mixing in a wide region and are barely relativis- 
tic by the end of the simulations. Models D2.5 and DOS 
have mixed deeply (the jet mass fraction is less than one 
everywhere) but keep larger Lorentz factors. Moreover, 
these models seem to have stopped the widening process 
of the mixing layer as it is deduced from the flattening 
of the mixing layer width as a function of time in Fig. [71 
Models BIO and DIO are also undergoing turbulent mix- 
ing. From Figs.dandjHI it is deduced that BIO and DIO 
are still mixing and transferring momentum by the end 
of simulations. These models will eventually lose a large 



amount of their initial longitudinal momentum, thereby 
becoming colder and denser due to mass cntrainmcnt from 
the ambient medium. Finally, models AlO, B20 and D20 
present a fast core ^ 1 Rj wide with rich internal structure 
as a consequence of the resonant modes (see subsection on 
the linear regime), which is surrounded by a hot and slow 
shear layer that extends up to 2 Rj in models AlO and 
B20 or 4 Rj in model D20. Let us point out that model 
AlO fFig. 115(1 displays a highly asymmetric structure, re- 
sulting from the development of resonant modes only on 
the upper interface. This is a consequence of the combina- 
tion of symmetric and antisymmetric modes, and probably 
of nonlinear interactions between resonant modes, which 
result in destructive interference on one side of the jet and 
constructive interference on the other. 
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Fig. 10. Snapshot at the last frame of the simulation of logarithmic maps of pressure, jet mass fraction and specific 
internal energy and non- logarithmic Lorentz factor for model A2.5. 
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Fig. 11. Same as Fig.Cnifor model B2.5. 



4. Discussion 

4.1. Non-linear stability 



ing models, confirmed the general trends resulting from 
the linear stability analysis: the faster (larger Lorentz fac- 
tor) and colder jets have smaller growth rates in the linear 



Simulations presented in Papers I and II, performed for regime. In Paper II, the non-finear evolution of the insta- 
the most unstable first reflection mode of the correspond- bility in the different models was characterized through 
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Fig. 12. Same as Fig.HUlfor model D2.5. 
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the processes of jet/ambient mixing and momentum trans- by the end of the simulation. Models in Class II (hot 

fer. The models were then classified into four classes (I to and fast jets) were slowly mixed in the nonlinear regime, 

IV) according to the particular nature of these processes progressively losing their longitudinal momentum. Models 

in each of the models. Class I models (corresponding to in Class III (hot and slow jets) have properties between 

cold and slow jets) were deeply mixed and mass-loaded Classes I and II. Finally, Class IV (containing cold/warm 



14 



M. Perucho et al.: Nonlinear stability of relativistic sheared planar jets 



Time = 500 R/c 



10 - 

6 - 



-6 - 
-10 ■ 



Log Pressure 



* I 4 « 
♦ Iff 



T.4BS>flf 



3.7ae-1B 



I 



Tracer 




5 io 15 so 25 ao 



5 10 IS so 25 30 



3.6Hi+01 



7.0S6+OB 



\ 



Log Specific Intamat Energy 
io\ 

6 - 



Raw Lorerrtz Factor 




S.B9ti-\-aO 



i.oa6+ao 




5 10 15 so 25 30 



5 10 IS so 25 30 
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Fig. 9. Evolution of the normalized total transversal mo- 
mentum in the jet as a function of time for all the simu- 
lations. Lines represent the same models as in Fig.El 

and fast models) appeared as the most stable in the non- 
linear regime. Shear layers formed in all the models as a 
result of the non-linear evolution. Models in Classes I/II 
developed broad shear layers and appeared totally mixed, 
cooled, and slowed down. In contrast, models in Classes 
III/IV have an inner core surrounded by thinner layers and 
keep a larger amount of their initial longitudinal momen- 
tum. We performed a number of additional simulations 
keeping the properties of the ambient medium fixed and 
changing the rest-mass density of the jet and the Lorentz 



factor. Results confirmed that these models behave like 
previous simulations, and are naturally placed in the clas- 
sification already defined. 

The stability classes considered in Paper II were de- 
fined according to the jet response to single modes. In this 
paper we revisit this classification scheme in the light of 
the present results based on more general perturbations. 
From the analysis of Figs. |H1 and El we classified jets 
depending on their nonlinear behaviour in three different 
groups: 

— Unstable 1 (USTl) models: jets which are disrupted 
after a strong shock is formed after the linear regime, 
enhancing turbulent mixing with ambient medium. It 
includes models A2.5, B2.5, D2.5, BOS, and DOS, 
i.e., lower Lorentz factor jets. The mixing layer width 
becomes larger than 5 Rj (Fig. [71), and they share more 
than 50% of their initial momentum with the ambient 
medium (Fig. (HI). 

— Unstable 2 (UST2) models: jets which are disrupted 
in the non-linear phase by a continuous process of mo- 
mentum transfer to the external medium, like BIO and 
DIO. This is observed in Fig. as a non-decreasing 
transversal momentum in the nonlinear regime. These 
models eventually end up sharing a large fraction of 
initial momentum and developing a wide mixing layer. 

— Stable (ST): jets which develop resonant modes and 
remain coUimated for long time, AlO, B20, and D20. 
These models have a thin mixing layer and share a very 
small fraction of their axial momentum with the am- 
bient medium. They expand, but remain coUimated. 
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In the course of their evolution, the jets develop a rich 
transversal structure in all the physical variables. This 
structure is different depending on the non-linear evolu- 
tion of the jets. Figure [201 displays the transversal profiles 
of relevant physical quantities averaged along the jet at 
the end of the simulations for model A2.5, representative 
of models in USTl, DIO of UST2, and B20 of ST. 

Model A2.5 shows a totally mixed, mass-loaded flow 
with averaged maximum speed 0.4c, i.e., barely relativis- 
tic, as these jets are efficiently slowed down by mass en- 
trainment after the disruption. The mass loading is in- 
ferred from the low values in the tracer profile (/ < 0.3), 
which imply a large fraction of ambient medium material 
inside the jet. The efficient conversion of kinetic energy 
into internal energy enhanced by the shock triggered in 
the early post-linear phase causes the jet to increase its 
specific internal energy. 

UST2 jets undergo a slower process of mixing, so they 
still keep a larger fraction of axial velocity and Lorentz fac- 
tor by the end of the simulation, even though they appear 
to be totally mixed (/ < 0.7 everywhere). However, as we 
have mentioned in the previous section, the mixing and 
slowing process is still going on in BIO and DIO, so it is 
clear that if the simulation had continued, the longitudinal 
velocity and Lorentz factor values would be smaller than 
those found. We also observe that the more mass-loaded 
parts of the jet (i.e., the region with — 10i?j < x < 0) are 
consistently colder. 

Finally, the jet in model B20 remains very thin. The 
velocity profile of the model has widened by 2 — 3 Rj by 



the end of the simulation, coinciding with the generation 
of a hot shear layer. This layer is seen in the figure as 
an overheated and underdense {p < 0.1) region shielding 
the unmixed core (/ = 1.), which keeps almost all its 
initial axial momentum and Lorentz factor. The core has 
a rich internal structure (see the pressure panel in Fig. I18|l 
that also manifests in the spiky structure of the shows 
longitudinal momentum profile. 

A comparison of the present non-linear evolution clas- 
sification scheme and that of Paper II (classes I-IV) allows 
us to conclude that, in general, models in Classes I and 
III fall into USTl, whereas models in Class II corresponds 
to UST2 and those in Class IV to ST. The reason models 
D2.5 and DOS (belonging formerly to Class III) move to 
USTl may be the inclusion of longer wavelength pertur- 
bations, along with the antisymmetric modes, which are 
more disruptive than the symmetric first reflection mode 
used in the previous work, and the lack of axial resolu- 
tion in the latter, as discussed in the Appendix of Paper 
II. This can be seen by comparing structures and evolu- 
tion of model D05 here and in Paper II, in particular the 
evolution of the mixing and momentum transfer. 

Regarding UST2 here compared to former class II, 
BIO and DIO undergo a very similar slow process of mo- 
mentum transfer to the external medium to that observed 
for DIO and D20 in Paper II, although their temperatures 
are very different and the shock in BIO is much stronger 
than in DIO (see Table 0)). The reason for this slow mo- 
mentum exchange may be the same as proposed in Paper 
II for models DIG and D20, i.e., a continuous conversion of 
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Fig. 16. Same as Fig.HUlfor model BIO. 

kinetic into internal energy due to the large initial Lorentz 
factor, which acts as a source of transversal momentum fa- 
voring the process of mixing and mass- loading. Model BIO 
changes from Class I in Paper II to UST2 here, meaning 
that disruption occurs by slow mixing in the new simula- 
tion, compared to sudden disruption in the previous one. 

Models in Class IV were characterized by a rich 
internal-structure jet preserving a large fraction of initial 
momentum and Lorentz factor. ST models share these fea- 
tures, but now we are able to clearly associate them to 
with the growth of resonant modes, which could be the 
reason for the breaking of the linear slope in model C20 
in Paper I (see Fig. 2 there). Steepening of short wave- 
length perturbations at the shear layer generates small 
shocks which favor local mixing and an efficient conver- 
sion of kinetic into internal energy. As a result of this 
process, the shear layer heats up and the jet expands to 
form a hot and underdense layer around the jet core (see 
Fig. l2()|l . It is remarkable that model AlO is largely asym- 
metric by the end of the simulation (see Fig. I15|) . This is 
a consequence of the resonant modes only growing on one 
side of the jet during the linear regime, and it is under- 
stood on the basis of asymmetry resulting from mixture of 
symmetric and antisymmetric modes. This effect, though 
much less evident, is also observed in model B20. Finally, 
model D20 has moved from class II in Paper II to ST here, 
clearly due to the appearance of resonant modes. This fact 
allows us to conclude that the fate of ST models would be 
exactly the same as those in UST2, if it were not for the 
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Fig. 21. Relativistic internal jet Mach number (Mj) ver- 
sus jet Lorentz factor (7) of the simulated models here and 
in Papers I and II. Symbols represent different non-linear 
behaviors: crosses stand for USTl disrupted jets (low rela- 
tivistic Mach number and low Lorentz factor); triangles for 
UST2 jets (moderately fast and supersonic), and squares 
for ST jets (highly supersonic and fast jets). Models with 
two different symbols are those with a different evolution 
in simulations presented here and those from Papers I and 
II (see text). 



growth of resonant modes; hence, their importance in the 
long term stability of these jets. 
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We classified jets depending on their nonlinear be- 
haviour in three different groups, which arc clearly sep- 
arated in a relativistic internal jet Mach number vs. jet 
Lorentz factor plot fFig. I21|l . In this plot, we also include 
models from Paper II, in order to show the general char- 
acter of our results and to note that this division of the 
stability properties of jets is more accurate than in Paper 
II, with the jet-to-ambient enthalpy ratio instead of the 
relativistic Mach number. A clear correlation between the 
two plotted parameters and the non-linear stability prop- 
erties of the jets is observed. Models BIO and D20 are 
represented with two different symbols to show the change 
of nonlinear behaviour from Paper II. These are placed 
in transition regions of the plot, either in Lorentz factor 
(BIO) or in relativistic Mach number (D20). This fact 
could explain differences in the non-linear behaviour as 
caused by changes in the initial jet profiles, what is quite 
evident in the case of D20, for resonant modes appear 
due to the presence of the shear layer. As in the previ- 
ous discussion, we have given the same symbols (crosses) 
as for USTl jets here to models in Class III of Paper II, 
as we do not consider that they have different non-linear 
behaviour in both simulations. Figure 1211 can be consid- 
ered as the relativistic counterpart of the M — ly (Mach 
number-density ratio) diagram in Bodo et al. (1994); note 
that the density ratio, v — pa/pj, is inverted with respect 
to ours. In our case the Mach number is relativistic; and 
the density ratio, which stands for the inertia of the flow, is 
replaced by the Lorentz factor here, as relativistic momen- 
tum is oc 7^, so it dominates the inertia of relativistic jets. 
Our conclusions are similar to theirs, for denser (higher 
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Lorentz factor) and highly supersonic jets (high relativis- 
tic Mach number) are the stablest. However, in our case, 
we found a higher degree of stability due to the growth of 
resonant, stabilizing modes. 

This result agrees with the conclusion of Hardee 
(2000), where hnear stability arguments show that dis- 
tortions induced by instabilities are smaller for higher 
Lorentz factor flows, although they were not associated 
to the shear resonances reported by us. 

Finally, simulations discussed in Appendices A (single 
antisymmetric mode in planar geometry) and B (single 
symmetric mode in cylindrical geometry) have confirmed 
general trends of the present clasification scheme, gener- 
alizing our results. 

4.2. Astrophysical implications 

One of the current open problems in extragalactic jet re- 
search is to understand the morphological dichotomy be- 
tween FRI and FRII jets. Several possible explanations 
have been proposed which mainly fall in one of these two 
general possibilities: either (i) FRI and FRII sources are 
intrinsically the same, and the morphology and jet evolu- 
tion depend mainly on the ISM in which they are embed- 
ded in the first kiloparsecs, or (ii) they depend on intrinsic 
differences stemming from the jet formation process (black 
hole rotation, Blandford 1994; accretion rate, Baum et al. 
1995; black hole mass, Ghisellini & Celotti 2001; the so- 
called magnetic switch, Meier et al. 1998), or (iii) a com- 
bination of both (e.g., Snellen & Best 2003). Of course, all 



18 



M. Perucho et al.: Nonlinear stability of relativistic sheared planar jets 



Time = 720 R/c 



Log Pressure 



Tracer 



3.046-02 





5 so 15 so SS 30 



5 10 15 so SS 30 




Fig. 18. Same as Fig.[Tnifor model B20. 



the mechanisms could come into play with differing effects 
and significance depending on the source. 

Leaving the basis of the morphological dichotomy 
aside, current models (Laing & Bridle 2002a,b and ref- 
erences therein) interpret FRI morphologies as the result 
of a smooth deceleration from relativistic (7 < 3, Pearson 
1996) to non-relativistic transonic speeds (~ 0.1 c) on kpc 
scales. On the contrary, flux asymmetries between jets 
and counter-jets in the most powerful radio galaxies and 
quasars indicate that relativistic motion (7^2 — 10) ex- 
tends up to kpc scales in these sources, although with 
smaller values of the overall bulk speeds (7 ~ 2 — 4, Bridle 
et al. 1994). Current models for high energy emission from 
powerful jets at kpc scales (e.g., Celotti et al. 2001) offer 
additional support to the hypothesis of relativistic bulk 
speeds on these scales. 

The results concerning the long-term evolution of rel- 
ativistic jets presented in this paper and summarized in 
Fig. 1201 confirm that slower and smaller Mach number 
jets (USTl) are entrained by ambient material and slowed 
down to w < 0.5 c after becoming overpressured (due to 
conversion of kinetic into internal energy) and being dis- 
rupted by nonlinear instabilities effects which cause flar- 
ing and rapid expansion of the mixing flow. UST2 jets 
undergo a smooth slowing down; and though by the end 
of the simulation jet velocity is ~ 0.9 c, this process is con- 
tinuous, and eventual loss of velocity to mildly relativistic 
values is to be expected. Finally, ST jets keep their ini- 
tial highly relativistic velocities, and their steadiness by 
the end of simulations makes them firm candidates for re- 



maining coUimated over long distances. Hence our results 
would point to a high Lorentz factor, highly supersonic 
jets as forming FRII Class, whereas FRI jets would be 
found in the opposite corner of the diagram (slow, small 
Mach number jets). The validity of our results extends 
to models with different jet-to-ambient-density ratios and 
specific internal energies as seen in Paper II. 

Our conclusions point to an important contribution 
by intrinsic properties of the source to the morphological 
dichotomy. Nevertheless, the importance of the ambient 
medium cannot be ruled out on the basis of our simula- 
tions, since we consider an infinite jet in pressure equilib- 
rium flowing in an already open funnel and surrounded 
by a homogeneous ambient medium. Thus our approach 
does not take into account the consequences of the in- 
teraction of the jet with the ambient in order to pene- 
trate it or the effects of a spatially varying atmosphere. 
Simulations following the spatial approach (perturbations 
grow with distance) for jets propagating in different ISM 
profiles and using a more realistic microphysics (allow- 
ing for a local mixture of electron, positron, and proton 
Boltzmann gases) will be performed in order to clarify 
these points. 

As dicussed in the introduction of this paper, there are 
plenty of arguments indicating the existence of transversal 
structure in extragalactic jets at all scales. In the simula- 
tions presented here, the initial states were defined with a 
continuous transition layer of thickness « Q.2Rj. As dis- 
cussed in the previous paragraphs, this shear layer has 
played an important role in the long-term stability of the 
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Fig. 19. Same as Fig.Hnifor model D20. 



jet flow. Besides this, thicker shear layers have been gener- 
ated in the course of the non-linear evolution. Relatively 
thin (« 2i?j), hot shear layers are found in present ST 
models (the physically meaningful counterparts of the lay- 
ers found in the three-dimensional, low-resolution simu- 
lations of Aloy et al. 1999, 2000), which could explain 
several observational trends in powerful jets at both par- 
sec and kiloparsec scales (see Aloy et al. 2000 and ref- 
erences therein). Conversely and according to our sim- 
ulations, these transition layers could be responsible for 
the stability of fast, highly supersonic jets, preventing the 
mass-loading and subsequent disruption. Finally, the type 
of shear layers developed by models USTl/2 could mimic 
the transition layers invoked in models of FRIs (Laing & 
Bridle 2002a,b). 

5. Conclusions 

We performed a number of simulations spanning a wide 
range of parameters, such as Lorentz factor and specific 
internal energy, for a general setup where a slab-sheared 
jet is perturbed with a set of symmetric and antisymmet- 
ric sinusoidal perturbations, in order to characterize the 
stability properties of relativistic jets. 

The most remarkable feature regarding the linear evo- 
lution of instabilities is the finding of resonant modes in 
our simulations, which were later confirmed by apply- 
ing linear stability theory to sheared flows. These modes 
are important for the long-term stability properties of 



some jets (ST), for they remain coUimated and unmixed, 
thereby keeping a large amount of initial axial momen- 
tum. Jets in which these modes do not grow fast enough 
with respect to longer modes are disrupted either after a 
shock or by slow momentum transfer and mixing. 

We classified jets depending on their nonlinear be- 
haviour in three different classes, which are clearly sep- 
arated in a relativistic internal Mach number vs. Lorentz 
factor plot (Fig. 1^ . USTl models are disrupted after 
a shock forms in the early post-linear phase, and ambi- 
ent gas penetrates deep into the jet stream, decelerating 
and cooling the initial flow down. UST2 models are slowly 
decelerated by an efficient conversion of kinetic energy 
into internal energy, which causes momentum transfer and 
mixing. Finally, ST models present little expansion, but 
remain collimated and isolated from the ambient by a hot 
shear layer. ST models would fall into UST2, if resonant 
modes were not present, as occurs for model D20 in Paper 
II. 

Our simulations admit a clear interpretation in the 
context of the morphological dichotomy of radio jets. 
Results presented here could point to high Lorentz factor, 
highly supersonic jets as forming FRII Class, whereas FRI 
jets would be related to slow, small Mach number jets. In 
the former, the generation of a hot shear layer surrounding 
a stable core could be related to the transversal structure 
observed in several powerful jets. 
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Fig. 20. Averaged transversal structure in the final state of the jets corresponding to models A2.5 (upper panels), 
DIO (middle), and B20 (bottom). Left panels (thermodynamical quantities): solid line, tracer; dotted line, rest mass 
density; dashed line, specific internal energy. Right panels (dynamical quantities): solid line, longitudinal velocity; 
dotted line, lorentz factor normalized to the initial value in the jet; dashed line, longitudinal momentum normalized 
to the initial value in the jet. Specific internal energy for model DIO was divided by 100 to fit in the scale. 
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